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Abstract 

Phase transitions of fluid mixtures of the type introduced by Stillinger and 
Helfand are studied using a continuum version of the invaded cluster algo- 
rithm. Particles of the same species do not interact, but particles of different 
types interact with each other via a repulsive potential. Examples of interac- 
tions include the Gaussian molecule potential and a repulsive step potential. 
Accurate values of the critical density, fugacity and magnetic exponent are 
found in two and three dimensions for the two-species model. The effect of 
varying the number of species and of introducing quenched impurities is also 
investigated. In all the cases studied, mixtures of g-species are found to have 
properties similar to g-state Potts models. 
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I. INTRODUCTION 



Several years ago Stillinger and Helfand introduced a simple but nontrivial model 
of fluid demixing. Their original model consists of a binary mixture of A and B particles. 
Particles of the same type do not interact with one another, but A and B particles interact 
with a repulsive potential such that the Mayer /-function is a Gaussian. This choice for the 
AB potential, known as the Gaussian molecule potential, greatly simplifies the calculation of 
virial coefficients and most work for this potential has been done using series methods |5],|J. 
The main motivation for this work is to confirm Ising universality for the critical exponents 
of continuum systems. 

In this paper we study the Stillinger- Helfand model and some of its generalizations using 
cluster Monte Carlo methods. Where possible, we compare our results to the series analyses 
and to results for the Ising-Potts universality classes. Although the Gaussian molecule 
potential yields a more tractable virial expansion, it is easier to implement a cluster algorithm 
for the repulsive step potential. We also consider a generalization of the Stillinger-Helfand 
model to q species (components), such that particles of the same species do not interact 
but particles of different species interact with a repulsive potential. We expect that this 
generalization will be in the same universality class as the g-state Potts model for q not too 
large and another motivation for this study is to confirm this correspondence. For example, 
we know that the two-dimensional (2D) Potts model for q > 4 has a first-order transition. 
Does the g-component 2D Stillinger-Helfand model also have a first-order transition for 
q > 4? In addition, we consider the effect of quenched disorder by randomly adding fixed 
scattering centers. There are general arguments || that quenched disorder causes first-order 
transitions to become continuous. These arguments hold rigorously for 2D Potts models || , 
but have not been studied for continuum models. 

In previous work cluster Monte Carlo methods were applied to the Widom-Rowlinson 
model || . The Widom-Rowlinson and Stillinger-Helfand models are closely related, the only 
difference is that the Widom-Rowlinson model has a hard-core interaction between different 
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species. In this paper, the invaded cluster Monte Carlo method introduced in Ref. is 
extended to soft-core repulsive potentials and is used to find the phase transition point for 
a given temperature without prior knowledge of the critical fugacity. The invaded cluster 
method has almost no critical slowing for the Widom-Rowlinson model, and we find that 
similar results hold for the Stillinger-Helfand models studied here. 

II. DESCRIPTION OF THE MODELS AND NOTATION 

We consider g-component (q > 1) fluids in d dimensions with d = 2, 3. The components 
(species) have no self-interaction but particles of one species interact with particles of all 
other species via an isotropic repulsive potential, U(r). We consider two choices for U(r), 

(U , if r < a 

Ktep(r) = ^ (2.1a) 
[ 0, if r > a , 

U gm (r) = -kT ln(l - e- r2/a2 ). (2.1b) 

The limit (3 = 1/kT — > oo for the step potential corresponds to the Widom-Rowlinson 
model. For the Gaussian molecule potential, the temperature T plays no role because the 
Boltzmann factor, e -1311 ^^ = 1 — e~ r ' 2 ^ a ' 2 is, by design, independent of T. 

In general, each component of the g-component fluid may have a distinct fugacity; how- 
ever symmetry considerations dictate that a demixing transition occurs with all fugacities 
equal, and hence we set all the fugacities equal to a single value, z. For sufficiently small 
z, the q species are mixed, while for large z, there are q distinct phases because different 
species repel one another, with each phase predominately composed of one species. If q is 
not too large, there is expected to be a single demixing transition separating these regimes 
that is in the same universality class as the g-state Potts model. 

For very large q, the correspondence between Potts models and Widom-Rowlinson mod- 
els must break down. Although Potts models have a single ordering transition, Widom- 
Rowlinson models can be presumed to have an intermediate crystalline phase for d > 3 and 
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large q. To understand this phase, consider the limits z< 1 and q ^> 1, with the product 
X = qz order unity. Then non-overlapping particles appear with an effective fugacity of 
A. However when two particles overlap, the cost is an additional factor of 1/q because the 
overlapping particles must be of the same species. Hence the limiting model is precisely 
the hard sphere gas which we presume has a crystalline phase in d > 3 0. It is therefore 
reasonable to assume that such a phase occurs in the Widom-Rowlinson models for large q 
and z of order g -1 . Needless to say, for fixed q, when the fugacity is sufficiently large, the 
model will demix, and thus the crystalline phase is an intermediate phase. We also expect 
an intermediate crystalline phase for large q soft-core Stillinger-Helfand models based on a 
mapping to a single component fluid with a repulsive soft-core potential. For example, the 



Gaussian core model |T0| is known to crystallize in d = 3. Although intermediate phases do 
not occur for the usual Potts models, they are not a consequence of the continuum; indeed 
such phases are known to occur on the lattice for the site-dilute (annealed) Potts models 



ITJ as well as for the lattice version of the Widom-Rowlinson model [T^ 



III. CLUSTER ALGORITHM 

The algorithms used here are, broadly speaking, examples of cluster algorithms of the 
type first introduced by Swendsen and Wang [|T3"l , |l4| . Cluster algorithms have been found to 
be much more efficient than local algorithms such as the Metropolis algorithm for simulating 
spin systems and lattice gases near critical points. Cluster algorithms would be very useful 
for off-lattice systems, but no general cluster method has yet been developed; indeed, the 
only off-lattice models for which highly efficient cluster methods are known are models of 
the Stillinger-Helfand and Widom-Rowlinson type. The distinguishing features of this class 
of models are that particles of the same species have no self-interaction and that there is 
a purely repulsive interaction between particles of different species. In this case, graphical 



representations and clusters algorithms are available fl5,|IE| and have been implemented for 
the Widom-Rowlinson model [TTfl. 
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Cluster algorithms for spin systems work by identifying clusters of spins and then ran- 
domly flipping these clusters. Cluster are defined by placing bonds between nearest neighbor 
aligned spins with a probability that depends on the temperature. For fluid systems, bonds 
are placed between particles of the same type with a probability that depends on the tem- 
perature and the interaction potential. Instead of flipping spins, clusters of particles are 
removed from the system and then new particles are added via a nonuniform Poisson pro- 
cess that depends on the fugacity and the potential due to the remaining particles. 

Cluster algorithms are typically used with fixed values of the external parameters such as 
temperature or fugacity. However, when the location of the phase transition is not known, 
much computational effort in studying the transition is spent in locating the transition. To 



avoid this problem, invaded cluster methods can be used which automatically adjust 

a thermodynamic parameter (for example, temperature or fugacity) to its value at the phase 
transition. This adjustment is accomplished by using the fact (proved for the q = 2 case 
fT5|| ) that the clusters just percolate at the transition. In invaded cluster algorithms, clusters 
are grown until a signature of percolation is observed. The value of the thermodynamic 
parameter at the transition is an output of the simulation obtained from the fraction of 
successful attempts to add particles or bonds to the system. The invaded cluster algorithm 
also may be used to distinguish first-order from continuous transitions as discussed in Ref. 
|T8| for Potts models. This method for distinguishing the order of the transition is discussed 
below and will be used in Section IIVC. 



We first describe the cluster algorithm discussed in Section 3.5 of Ref. JT5[ for Stillinger- 
Helfand models and then discuss how it can be modified to be an invaded cluster algorithm. 
We assume that we have a configuration consisting of particle positions and a set of bonds 
connecting some of the particles and describe how to obtain the next configuration: 

1. Identify all clusters of particles defined by the bonds. A particle with no bonds is con- 
sidered to be a singleton cluster. For each cluster, independently and with probability 
1/q, label it a black cluster and with probability 1 — 1/q label it white. 
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2. Remove all particles in black clusters. The remaining white particles are at a set of 
positions W. 

3. Replenish the black particles via a Poisson process with local intensity y(x) given by 

y {x) = ze-W x \ (3.1a) 

V(x)= J2U(\x-y\). (3.1b) 

yew 

where z is the fugacity and U (r) the potential. 

4. For each pair of black particles, place a new bond between them with probability p(r) 
given by 

p(r) = l-e- mr \ (3.2) 

where r is the separation between the particles. Note that p(r) is minus the Mayer 
/-function for the potential. 

5. Eliminate the white and black labels for the clusters. 

This procedure comprises one Monte Carlo step. 

Given a configuration of particle positions and bonds without species labels, it is possible 
to obtain a full multicomponent configuration where each particle has a species label. This 
assignment is accomplished by identifying clusters and then randomly and independently 
assigning one of the q species labels to each cluster. The species label of each particle is 
taken to be the species label of its cluster. This labeling of particles is only possible if q 
is a positive integer. However, the algorithm makes sense for all q > 1, in analogy to the 
relation between Potts models, which are defined for positive integer q, and random cluster 
models which interpolate between them and are defined for all q > 1. 

It is instructive to consider the nature of the cluster configurations generated by the 
algorithm as a function of the fugacity for a fixed temperature. Suppose that the fugacity 
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and, hence, the density is very small. Then p(r) is typically small because the particles 
are far apart, and most clusters are singletons. In step 2, a fraction 1/q of the particles is 
removed. In step 3, particles are replenished as a nearly ideal gas because the exponential 
factor in Eq. (|3.1a| ) is a small perturbation except in the vicinity of the remaining particles. 
The end result is a nearly ideal multi-component gas. In the limit of large fugacity and 
density, we expect a phase in which a single species is predominant with a small admixture 
of the other species. The bonds connecting particles of the dominant species are sufficiently 
dense that almost all members of this species are in a single large cluster. The minority 
species are almost all in widely scattered singleton clusters. When the majority species is 
white, as occurs in about 1/q of the Monte Carlo steps, the large cluster is removed and 
then replaced as a nearly ideal gas in a slightly perturbed background potential generated 
by the minority species. An important feature of this picture is that the clusters do not 
percolate at small fugacity and do percolate at large fugacity. At some intermediate value 



of the fugacity, there must be a percolation transition. As discussed in Refs. |Lo[[L9[gO] , the 
percolation transition of the clusters coincides with the demixing transition of the fluid. 

The coincidence of the percolation transition and the demixing transition justifies an 
invaded cluster version of the above cluster algorithm. The invaded cluster algorithm is very 
similar to the fixed z cluster algorithm described above except that steps 3 and 4 are modified 
as follows. Instead of putting down new black particles as a Poisson process at a fixed 
intensity, black particles are added to the system one at a time according to the potential 
V(x) (see below). After each black particle is added, bonds between the new particle and 
all previously placed black particles are put down with probability p(r) = 1 — exp(— /3U(r)). 
The black clusters defined by these bonds are monitored after each particle is added, and the 
process of adding particles is stopped when a stopping condition is satisfied. For simulating 
the phase transition, the stopping condition is that one cluster spans the system. For periodic 
boundary conditions, spanning is taken to mean that a cluster wraps around the system in 
at least one of the d directions. The spanning condition insures that the algorithm simulates 



the phase transition [18 



In practice, a particle is added to the system according to the potential V(x) by the 
following procedure. A particle is tentatively placed at a random position x. A random 
number r is chosen in the interval [0, 1), and the particle placement is accepted if 

r < e~ mx) ; (3.3) 

otherwise the particle is rejected and another attempt is made to place a particle. Let 
z = (Ntot/L ), where iVtot is the total number of attempted particle placements in a Monte 
Carlo step, including both accepted and rejected placements, L d is the system volume, and 
the brackets (...) indicate an average over the simulation. Because the intensity y(x), defined 
in Eq. (|3.1a|) , and the Boltzmann factor e~^ vl ^) governing particle placements differ by a 
factor of the fugacity, we conclude that z is an estimator of z c , the value of the fugacity at 
the transition. Note that if the fluctuations cr§ in N tot /L d are small, then the invaded cluster 
algorithm is essentially identical to the fixed fugacity algorithm operating at z — z c . This 
identification justifies the use of the invaded cluster method. A more complete discussion of 
the invaded cluster method and the use of z as an estimator of a critical parameter is given 
in Ref. M. 



Whenever the invaded cluster method simulates a system at its critical point, scaling 
methods can be used to obtain critical exponents from the size dependence of divergent 
thermodynamic quantities such as the compressibility or the susceptibility. To study the 
latter, we consider the quantity 

X s p<E*?> (3-4) 
where Sj is the number of particles in the zth cluster. We now show that x is related to the 
usual susceptibility. Consider, for simplicity, the discretized version of the Stillinger-Helfand 
model on a lattice of linear dimension L with spacing e so that the total number of sites is 
[L/e] d . The demixing order parameter at site x is given by Spi(x) = ri\{x) — n(x)/q, where 
rii(x) = 1 if there is a particle of type 1 at site x and n\(x) = otherwise; n(x) counts the 
presence of a particle of any type. The relevant susceptibility x is defined by the second 
derivative of the pressure with respect to the (ordering) chemical potential: 
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X=- R J2(S P i(x)5 Pl (y)). (3.5) 

L x,y 

(The reason that e does not enter explicitly into Eq. ( |3.5| ) is that the derivatives are with 
respect to the log of the activity and it is the activity that is scaled by e.) For a given 
particle and bond configuration, averaging over assignments of species labels, it is clear that 
5pi(x)5pi(y) vanishes unless the sites x and y are both occupied and in the same cluster, in 
which case the result is q~ 2 (q — 1). Thus, for a fixed particle-bond configuration, we obtain 
the number of particles in the cluster at x if we sum over y. Summing over x yields the sum 
of the squares of the cluster sizes so that x — <?~ 2 (<? ~~ l)x> an d hence we conclude that x is 
related to the usual susceptibility. Finally, finite size scaling predicts that 

X ~ Vl\ (3.6) 

so that the scaling of x with system size can be used to extract the magnetic exponent 7. 

Cluster methods also may be used to distinguish first-order from continuous transitions. 
For this purpose, a fixed density stopping rule is used. Black particles are added to the 
system until the density p reaches a fixed value and then z is measured. In this way the 
canonical ensemble is simulated rather than the grand canonical ensemble. This procedure 
is done for a range of densities near the transition. If the transition is continuous, the 
fugacity is a strictly increasing function of p. However, if the transition is first-order, then 
the fugacity does not increase monotonically with increasing p in the coexistence region. 
Why does the nature of the z versus p curve signify whether a transition is continuous or 
first-order? Suppose that the demixing transition of a g-component system is first-order. At 
the transition, there is coexistence of q + 1 phases; q demixed phases and one mixed phase. 
Because the repulsive interaction is reduced for the demixed phases, these phases have a 
higher density than the mixed phase. Thus, in the thermodynamic limit, there is a range of 
p for which the fugacity is constant. Let pi be the density of the mixed phase and p2 the 
density of the demixed phase. Because In z = ds/dp, where s is the entropy density, we have 
that s is a linear function of p in the coexistence region. More specifically, s(p) is a linear 
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combination of s(pi) and s(p2), the entropy densities of the mixed and demixed phases. 
The linearity of s(p) applies in the thermodynamic limit. However, for a finite system, 
the entropy density is not linear in the coexistence region. Consider a system with linear 
dimension L and periodic boundary conditions at density p. This system also can be viewed 
as an infinite system with periodic constraints on the particles. Let s(p,L) be the entropy 
density of this periodically constrained system. Now suppose the constraints are removed 
and the system comes to equilibrium. If p\ < p < p 2 , demixing will occur spontaneously so 
that s(p, L) > s(p) with the equality holding only at the endpoints of the coexistence range. 
Because Inz = ds/dp, we must have that z is non-monotone in the coexistence region. This 
approach for distinguishing the order of a transition is very similar to the microcanonical 



Monte Carlo method used in Ref. 21 



IV. RESULTS 



In Section [IV A| we present results for the 2D and 3D Stillinger-Helfand Gaussian molecule 
models. The two-component step potential model is discussed in Section [IV B| and the q- 
component step potential is discussed in Section |1V C[ 



A. Gaussian molecule model in two and three dimensions 

We simulated the Gaussian molecule model (with the potential U gm defined in Eq. ( |2.1b| )) 
using the invaded cluster method and the spanning rule described in Section |T| for a range 
of linear dimensions L up to 140 in d = 2 and 40 in d = 3. We choose units such that 
distances are measured in units of a. We collected statistics for the number of particles 
in the spanning cluster M, the critical density p, the susceptibility x, the estimator of the 
critical fugacity z and its standard deviation <7g, and the normalized autocorrelation function 
for the spanning cluster size, Ym- For each value of L we averaged over 10 5 Monte Carlo 
steps. The estimator of the critical density is the average number of particles (of any species) 
per unit area (volume) when the spanning condition is fulfilled. Although U gm (r) does not 
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go to zero at finite r, it becomes very small for larger r and to speed the calculation, we set 
Ugm(r) = for r > 3. 

Tables | and y show the L dependence of M, p, \i %i a z, and tm for the 2D and 3D 
Stillinger-Helfand models, respectively. The integrated autocorrelation time tm is defined 
by 

1 oo 

^ = , + Er M (t). (4.i) 
z t=i 

This time is approximately the number of Monte Carlo steps between statistically indepen- 
dent configurations and enters into the error estimate for M. In practice, becomes 
indistinguishable from the noise for t ~ 10 Monte Carlo steps, and it is necessary to cut off 
the upper limit of the sum defining t m when the magnitude of T M becomes comparable to 
its error. 

Note that the fluctuations a z in z decrease with increasing L and that tm is small and 
hardly increases with L. These results demonstrate the validity and efficiency of the invaded 
cluster algorithm. The decrease in a z shows that as L increases, the invaded cluster becomes 
essentially equivalent to a fixed parameter cluster algorithm for which detailed balance can 
be proven. 

The error estimates for all quantities in Tables | and |T| except tm were obtained by 
computing the standard deviation of the quantity of interest and dividing by the square root 
of the number of measurements. This error estimate does not take into account correlations 
between successive Monte Carlo steps. To account for correlations, the error estimates in 
the tables for an observable O must be multiplied by \f2r , where To is the integrated 
autocorrelation time for O. The statistical errors for quantities derived from fits such as p c 
and j/u include the factor a/2t except that t q is replaced by t m . 

Figure [I] shows the results for p(L) versus 1/L for the 2D Stillinger-Helfand model. The 
value of p in the limit of L —>■ oo is extrapolated from the finite size data by doing a linear 
least squares fit omitting the values for L = 20 and 40 yielding the result, p c (2) = 1.1644 ± 
0.0004. A similar extrapolation for the critical fugacity yields z c (2) = 1.3536 ± 0.0008. 
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Similarly, extrapolating the result for the 3D Stillinger-Helfand model using the data for all 
available L yields p c (3) = 0.440 ± 0.001 and z c (3) = 0.5826 ± 0.0013. Our error values for 
these critical parameters are one standard deviation from the linear least squares fit of the 
fugacity or density versus 1/L; no effort has been made to estimate systematic errors. All 
of the fits have acceptable goodness-of-fit probability values Q. 

Our 3D value for the critical density is consistent with the series result of Lai and Fisher 
J| ; p c (3) = 0.441 ± 0.001 (Eq. (36) of Ref. 0) but our critical fugacity is somewhat larger 
than their value, z c (3) = 0.5785 ± 0.0002 (Eq. (44) of Ref. @). Note that Lai and Fisher 
report results using a different convention so that their values of p c and z c must be divided 
by ix d l 2 to compare with our values. 

The exponent ratio 7/V can be obtained from the scaling of the susceptibility x with 
L according to Eq. (|3.6|) . Figure |2| shows a log- log plot of x versus L for the 2D Gaussian 
molecule model. A least squares fit of all the data to a simple power law does not yield 
an acceptable goodness of fit value Q. If the smallest value of L is omitted, we obtain 
7/1/ = 1.745 ± 0.001 with x 2 = 6.2, Q = 0.19, and DF = 4 (degrees of freedom). The Q 
value indicates a reasonable fit to a simple power law, but the fitted value of 7/1/ is 5 a from 
the 2D Ising value of 7/V = 7/4. A reasonable explanation of this result is that the 2D 
Gaussian molecule is, indeed, in the 2D Ising universality class, but that there are relatively 
slowly varying corrections to scaling. 

A least squares fit to all the data for the susceptibility x f° r the 3D Gaussian molecule 
model yields 7/// = 1.9626 ± 0.0044 with x 2 = 0.11, DF = 2, and Q = 0.95. The Q value 
near unity suggests that the data is well fit to a pure power law. Recent high precision 
Monte Carlo studies of the 3D Ising model [^] yield 7/1/ = 1.9630(30) which is consistent 
with our results. Our results add weight to the hypothesis that the Stillinger-Helfand model 
is in the Ising universality class for both 2D and 3D. The relatively high precision results for 
7/V from the Gaussian molecule model suggests that models of the Stillinger-Helfand type 
may be useful for high precision studies of the 3D Ising universality class. The isotropy of 
the interaction and absence of an underlying lattice might make for smaller corrections to 
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scaling in Stillinger-Helfand models compared to lattice spin models. 



B. Step potential in 2D 

Table [ITJ summarizes our results for the 2D Stillinger-Helfand model with the step po- 
tential given by Eq. ( J2.1a ) and temperature T = 1 (measured in units of Uq). For each 



value of L we averaged over 10 Monte Carlo steps. The results are qualitatively similar to 



the Gaussian molecule model. Table [TV] shows the temperature dependence of the measured 
quantities for L fixed at L = 20. The values of p and z at low temperatures should reduce 
to the Widom-Rowlinson model. In Ref. we measured the critical parameters of the 
Widom-Rowlinson model using the invaded cluster method. For L = 40 (the smallest size 
measured) we obtained p = 1.525 and z = 1.720, values that are close to the values of p 
and z for the two lowest temperatures in Table |V]. This agreement confirms that the step 
potential is continuously connected to the hard core potential. 

If the L = 20 data point is omitted, we obtain from a least squares fit to the data for 
the susceptibility x, ll» = 1-7434 ± 0.0009 with X 2 = 0.52, Q = 0.81, and DF = 3. 

C. Dependence of the order of the transition on q and on impurities 

The critical properties of the g-component Stillinger-Helfand model are expected to 
be closely related to the g-state Potts model. One of the features of the g-state Potts 
model is that the transition is continuous for small q and is first-order for q > q c {d), where 



q c (2) = 4 and 2 < g c (3) < 3. We have used the method described in Section |TTT| to determine 
the order of the transition as a function of q for the g-component Stillinger-Helfand step 
potential model. Figure |3| shows the fugacity z as a function of p for d = 2 for L = 40 and 
T — 1. Note that for q = 3 the curve is clearly monotonically increasing, which implies 
a continuous transition. For q > 5 the curves are clearly non-monotonic, which implies a 
first-order transition. For q = 4 the curve is essentially flat within the error bars (whose 
size is approximately that of the symbols). Although the effective value of g c is expected 
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to vary with L, these results are consistent with the hypothesis that q c {2) = 4 for the 2D 
Stillinger-Helfand step potential model. 

Figure ^ shows z as a function of p for the 3D Stillinger-Helfand step potential model 
for L = 20 and T — 1. The q = 2 curve is clearly monotonically increasing while the q = 3 
is clearly not, implying that 2 < q c (3) < 3 as for the 3D Potts model. 

Finally, we have studied the effect of quenched impurities on the nature of the transition 
for the q = 3 Stillinger-Helfand step potential model. The impurities consist of randomly 
placed scatterers that interact with all the fluid particles via the same repulsive step potential 
that exists between different components. Figure |5| shows a plot of z versus p for the q = 3 
Stillinger-Helfand step potential model in 3D for four impurity densities ranging from 0.025 
to 0.0625. For each of the 10 impurity configurations considered for a given density, data from 
10 3 Monte Carlo steps are collected. For the two lowest impurity concentrations, the z versus 
p curve is non-monotonic as is the case for the pure system, while for the two highest impurity 
concentrations, the curve is monotonic indicating a crossover to a continuous transition. This 
behavior is in accord with general arguments that the presence of quenched impurities 
should cause a first-order transition to become continuous. It is not clear from our data 
whether there is a critical value of the disorder below which the transition remains continuous 
or whether the crossover at finite disorder strength is a finite size effect and that any strength 
of disorder is sufficient to make the transition continuous in the thermodynamic limit. 

V. DISCUSSION AND CONCLUSIONS 

We have studied the Stillinger-Helfand model and several generalizations using the 
invaded cluster algorithm. Our results for g-component Stillinger-Helfand models with 
2 < q < 8 are consistent with the hypothesis that these models are in the same universality 
class as the corresponding Potts models. In addition, we have shown that the addition of 
quenched disorder causes the demixing transition to change from first-order to continous 
for those values of q for which the pure system transition is first-order. For the case q = 2 
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and d — 2, our results for the magnetic exponent are outside the statistical error bars of 
the exact Ising value. However, we believe that this difference is most likely the result of 
slowly varying corrections to scaling. It would be useful to consider larger system to confirm 
Ising universality. It would also be interesting to consider larger values of q to explore the 
possibility of an intermediate crystalline phase in Stillinger-Helfand models. 



VI. ACKNOWLEDGEMENTS 

This work was supported by NSF grants PHY-9801878 (RS), DMR-9633385 (HG), DMR- 
9978233 (JM), DMS-9971016 (LC), and NSA grant MDA904-98-1-0518 (LC). We thank 
Gregory Johnson for useful discussions. 

^Present address: Courant Institute of Mathematical Sciences, 251 Mercer Street, New York, 
NY 10012. 



15 



REFERENCES 
[1] F. H. Stillinger and E. Helfand, J. Chem. Phys. 41, 2495 (1964). 
[2] E. Helfand and F. H. Stillinger, J. Chem. Phys. 49, 1232 (1968). 
[3] A. Baram, M. W. Maddox, and J. S. Rowlinson, Mol. Phys. 76, 1093 (1992). 
[4] S.-N. Lai and M. E. Fisher, Mol. Phys. 88, 1373 (1996). 

[5] K. Hui and A. N. Berker, Phys. Rev. Lett. 62, 2507 (1989); 63, 2433 (E) (1989). 
[6] M. Aizenman and J. Wehr, Phys. Rev. Lett. 62, 2503 (1989). 

[7] G. Johnson, H. Gould, J. Machta, and L. K. Chayes, Phys. Rev. Lett. 79, 2612 (1997). 

[8] B. Widom and J. S. Rowlinson, J. Chem. Phys. 52, 1670 (1970). 

[9] H. L. Frisch and J. K. Percus, Phys. Rev. E 60, 2942 (1999). 
[10] F. H. Stillinger and T. A. Weber, J. Chem. Phys. 68, 3837 (1978). 
[11] L. Chayes, R. Kotecky, and S. Shlosman, Comm. Math. Phys. 171, 203 (1995). 
[12] L. K. Runels and J. L. Lebowitz, J. Math. Phys. 15, 1712 (1974). 
[13] R. H. Swendsen and J. S. Wang, Phys. Rev. Lett. 58, 86 (1987). 

[14] M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in Statistical Physics, 
(Oxford U. Press, Oxford, 1999). 

[15] L. Chayes and J. Machta, Physics A 254, 477 (1998). 

[16] O. Haggstrom, M. N. M. Lieshout, and J. Moller, Bernoulli 5, 641 (1999). 

[17] J. Machta, Y. S. Choi, A. Lucke, T. Schweizer, and L. V. Chayes, Phys. Rev. Lett. 75, 
2792 (1995). 

[18] J. Machta, Y. S. Choi, A. Lucke, T. Schweizer, and L. M. Chayes, Phys. Rev. E 54, 
1332 (1996). 

16 



[19] W. Klein, Phys. Rev. B 26, 2677 (1982). 

[20] J. T. Chayes, L. Chayes, and R. Kotecky, Comm. Math. Phys. 172, 551 (1995). 

[21] D. H. E. Gross, A. Ecker, X. Z. Zhang, Annalen Phys. 5, 446 (1996). 

[22] H. W. J. Blote, E. Luijten, and J. R. Heringa, J. Phys. A: Math. Gen. 28, 6289 (1995). 



17 



TABLE I. The L dependence of the number of particles in the spanning cluster M, the critical 
density p, the susceptibility x, the estimator of the critical fugacity z, its standard deviation as, 
and decorrelation time tm for the 2D Gaussian molecule potential. The averages are over 10 5 
spanning clusters. The error estimates were obtained by computing the standard deviation of the 
quantity of interest and dividing by the square root of the number of measurements. The error 
estimates for tm are obtained from the variation of r with the upper limit in the summation of 
Eq. ( |4.lD and hence represent an estimate of the systematic error. The autocorrelation function 
r^f (i) is distinguishable from the noise for t ~ 10 Monte Carlo steps. 



L 


M 


P 


X 


z 






20 


297.1(1) 


1.105(1) 


245.2(5) 


1.3286(8) 


0.257 


0.56(1) 


40 


1095(1) 


1.1317(4) 


832.5(10) 


1.3469(6) 


0.183 


0.63(3) 


60 


2346(2) 


1.1418(3) 


1695(2) 


1.3504(5) 


0.149 


0.68(3) 


80 


4017(4) 


1.1472(3) 


2796(2) 


1.3516(5) 


0.131 


0.71(5) 


100 


6098(5) 


1.1503(2) 


4119(3) 


1.3519(4) 


0.117 


0.75(3) 


120 


8587(7) 


1.1532(2) 


5670(3) 


1.3524(4) 


0.108 


0.76(3) 


140 


11462(10) 


1.1547(2) 


7419(5) 


1.3519(4) 


0.0982 


0.79(2) 
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TABLE II. The L dependence of M, p, x, z, erg, and tm for the 3D Gaussian molecule potential. 
The averages are over 10 5 spanning clusters. The error estimates are calculated as discussed in 
the caption of Table |. represents one standard deviation. The autocorrelation function Tj^(t) is 
distinguishable from the noise for t ~ 60 Monte Carlo steps. 



L 


M 


P 


X 


z 






10 


171.6(3) 


0.435(2) 


35.7(2) 


0.577(2) 


0.103 


0.52(4) 


20 


956(2) 


0.438(1) 


138.9(4) 


0.580(1) 


0.0582 


0.55(3) 


30 


2614(3) 


0.438(1) 


308(1) 


0.581(1) 


0.0424 


0.57(2) 


40 


5344(6) 


0.439(1) 


542(2) 


0.581(1) 


0.0341 


0.58(3) 



TABLE III. The L dependence of M, 


P, X, z, a 


z and tm for the 2D Stilling 


jer-Helfand step 


potential model at T = 


: 1. The averages 


are over 


10 6 spanning clusters. 


The autocorrelation 


function Tm^) is distin^ 


mishable from the 


noise for t 


~ 10 Monte Carlo steps. 




L M 


P 


X 


z 






20 509.5(2) 


1.9249(3) 


720.7(5) 


2.2626(4) 


0.379 


0.585(5) 


40 1871.8(6) 


1.9633(2) 


2429(2) 


2.2819(3) 


0.270 


0.665(10) 


60 4001(2) 


1.9779(3) 


4927(3) 


2.2856(4) 


0.219 


0.702(10) 


80 6855(2) 


1.9855(2) 


8131(4) 


2.2865(3) 


0.191 


0.737(10) 


100 10411(3) 


1.9902(1) 


12000(6) 


2.2864(2) 


0.169 


0.772(15) 


120 14651(4) 


1.9934(1) 


16494(7) 


2.2860(2) 


0.154 


0.783(10) 
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TABLE IV. The T dependence of M, p, x, z, <?z, and tm for the 2D Stillinger-Helfand step 
potential model at L = 20. The averages are over 10 6 MC steps. The autocorrelation function 
rjtf(f) is distinguishable from the noise for t ~ 60 Monte Carlo steps. The computational time 
required for a L = 20 system of T = 10 is approximately a week on a 533 MHz Alpha processor. 



T 


M 


P 


X 


z 


<y~z 




0.05 


411.1(2) 


1.4905(3) 


470.3(3) 


1.7015(4) 


0.310 


0.58(0.5) 


0.2 


412.2(2) 


1.4950(2) 


472.4(3) 


1.7072(4) 


0.304 


0.57(0.5) 


0.5 


439.1(2) 


1.6076(2) 


535.9(3) 


1.8508(4) 


0.326 


0.58(0.5) 


1 


509.5(2) 


1.9249(2) 


720.6(4) 


2.2626(4) 


0.379 


0.59(1) 


3 


763.0(3) 


3.3065(4) 


1612.7(8) 


4.1052(6) 


0.598 


0.62(0.5) 


5 


963.1(7) 


4.6606(4) 


2567(2) 


5.9353(8) 


0.785 


0.64(1) 


7 


1133(1) 


5.9934(5) 


3551(2) 


7.7433(10) 


0.957 


0.65(1) 


10 


1353(1) 


7.9680(9) 


5059(4) 


10.426(2) 


1.182 


0.66(1) 
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0.01 0.02 0.03 0.04 0.05 

1/L 

FIG. 1. Plot of p versus 1/L for the 2D Stillinger-Helfand Gaussian molecule model. The 
straight line represents a least squares fit omitting the value for L = 20 and 40. The extrapolated 
result is p c = 1.1644 ± 0.0004. 
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ln(L) 

FIG. 2. Log-log plot of x versus L for the 2D Gaussian molecule model. The straight line 
represents the best fit to the data, omitting L = 20, with 7/1/ = 1.745. 
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FIG. 3. Plot of the fugacity z versus p for the 2D Stillinger-Helfand step potential model at 
T = 1, L = 40, for q = 3,4,5,6, and 8. Each point is averaged over 20,000 MC steps. The error 
bars are of the size of the markers. 
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FIG. 5. Plot of z versus p for the q = 3 Stillinger-Helfand step potential model in 3D with 
quenched impurities at T = 1 and L = 20. The 4 traces in the graph correspond to 200, 300, 400 
and 500 fixed impurity particles, corresponding to impurity densities equal to 0.025, 0.0375, 0.05, 
and 0.0625. 
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